function out = getExpectedExcessReturns(loadData,setup,output,theta2,h)

% We consider all maturities
tmp                = output.model;
setupAll           = setup;
setupAll.matSelect = 1:max(setup.matSelect);
if output.model.modelNum == 1
    [modelAll,errorMes]= solveQTSM(tmp.alpha,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setupAll);
elseif output.model.modelNum == 2
    [modelAll,errorMes] = solveQTSM_constRsmooth(tmp.alpha,tmp.rhor,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setupAll);
end
% Adding vecC to model
modelAll.vecC = zeros(length(modelAll.A),modelAll.nx^2);
for k=1:length(modelAll.A)
    modelAll.vecC(k,:) = reshape(modelAll.C(:,:,k),1,modelAll.nx^2);
end
modelAll.modelNum = output.model.modelNum;
modelAll.factorSignResOn = 0;
modelAll.macros   = output.macros;
modelAll.rLag     = output.model.rLag;
modelAll.rhor     = output.model.rhor;
econAct           = setup.econAct;
numSim            = 1D4;

% Realized excess holding period return in the data
maxMat           = max(modelAll.matSelect);
annualMaturities = modelAll.matSelect(1:maxMat)/12;
yields           = loadData.yields_h1(2:end,1+(1:maxMat)); 
if h == 1
    shortRate = [loadData.crspRiskFree_h1(2:end,1);nan];
elseif h == 3
    shortRate = [loadData.crspRiskFree_h3(2:end,1);nan];
end
[T,numYields]    = size(yields);
% Note that we loose the last h observations, because we cannot compute excess returns at T-h+1,T-h+2,...,T
holdPerReturn    = [nan(T-h,h),...
                    (repmat(annualMaturities(1,h+1:end),T-h,1).*yields(1:T-h,1+h:end)...
                    -repmat(annualMaturities(1,1:end-h),T-h,1).*yields(1+h:T,1:end-h))];
xhrData          = holdPerReturn/h - repmat(shortRate(1:T-h,1)/12,1,numYields);  
                
% For the P dynamics we use estimates from Step 3. Note that we
% here apply the sigma from step 3
factors = [output.macros'; output.factors];
nx      = size(factors,1);
[h0_1,h0_2,hx_1,hx_2,sigma,gama0,gamaz,gamax,sigmaz] = unfoldTheta2_threshold_macro(theta2,nx);

% % The shocks for forecasting the states
rng(1,'twister')
shocks        = randn(nx,numSim,h);
shocks        = shocks-repmat(mean(shocks,2),1,numSim,1);    %ensure conditional zero mean
shocks        = shocks./repmat(std(shocks,1,2),1,numSim,1);  %ensure conditional unit variance
shocksz       = randn(numSim,h);
shocksz       = shocksz - repmat(mean(shocksz,1),numSim,1);  %ensure conditional zero mean
shocksz       = shocksz./repmat(std(shocksz,1,1),numSim,1);  %ensure conditional unit variance
xhrModel      = zeros(T-h,numYields);
tmpR          = nan(numSim,h);
for t=1:T-h
    x                = factors(:,t);
    z                = econAct(t,1);
    factorSim        = nan(nx,numSim,h+1);
    factorSim(:,:,1) = repmat(x,1,numSim);
    zSim             = nan(numSim,h+1);
    zSim(:,1)        = z;
    for i =1:h
        boom  = (zSim(:,i) >= 0)';
        reces = (zSim(:,i) < 0)';
        factorSim(:,:,i+1) = repmat(boom,nx,1).*repmat(h0_1,1,numSim) ...
            + repmat(reces,nx,1).*repmat(h0_2,1,numSim) ...
            + repmat(boom,nx,1).*(hx_1*factorSim(:,:,i)) ...
            + repmat(reces,nx,1).*(hx_2*factorSim(:,:,i))...
            + sigma*shocks(:,:,i); %randn(nx,1);
        zSim(:,i+1) = repmat(gama0,numSim,1) + gamaz*zSim(:,i) ...
            + (gamax'*factorSim(:,:,i))' + sigmaz*shocksz(:,i); %randn(1,1);
        
        if nx == 3
            xxTmp = [factorSim(1,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(2,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(3,1:numSim,i+1).*factorSim(:,1:numSim,i+1)];
        elseif nx == 4
            xxTmp = [factorSim(1,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(2,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(3,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(4,1:numSim,i+1).*factorSim(:,1:numSim,i+1)];
        elseif nx == 5
            xxTmp = [factorSim(1,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(2,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(3,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(4,1:numSim,i+1).*factorSim(:,1:numSim,i+1); ...
                factorSim(5,1:numSim,i+1).*factorSim(:,1:numSim,i+1)];
        else
            error('in getExpectedExcessReturn, nx > 5');
        end
        tmpR(:,i) = modelAll.r0 + modelAll.rx*factorSim(:,:,1+i) + modelAll.rxx*xxTmp;
    end
    
    if modelAll.modelNum == 1
        ExpYields_h = mean(modelAll.g0 + modelAll.gx*factorSim(:,:,h+1) + modelAll.gxx*xxTmp,2);
        yields      = modelAll.g0 + modelAll.gx*x + modelAll.gxx*kron3(x,x);
    elseif modelAll.modelNum == 2       
        if h == 1
            ExpYields_h = mean(modelAll.g0 + modelAll.rhor*modelAll.rLag(t,1) + modelAll.gx*factorSim(:,:,h+1) + modelAll.gxx*xxTmp,2);
        else
            rExp(1,1) = modelAll.r0 + modelAll.rhor*modelAll.rLag(t,1)...
                       +modelAll.rx*x + modelAll.rxx*kron3(x,x);
            for i=2:h
                rExp(i,1) = modelAll.rhor*rExp(i-1,1)+ mean(tmpR(:,i),1);
            end
            % Note that rExp(h,1) is r_t+h-1 
            ExpYields_h = modelAll.g0 + modelAll.rhor*rExp(h,1) + mean(modelAll.gx*factorSim(:,:,h+1) + modelAll.gxx*xxTmp,2);
        end
        yields      = modelAll.g0 + modelAll.rhor*modelAll.rLag(t,1) ...
                    + modelAll.gx*x + modelAll.gxx*kron3(x,x);
    end
    % xhr_t,h^(n) = (n*y_t^(n)-(n-h)*E_t[y_t+h^(n-h)])/h-y_t^(h)/12
    tmp           = (annualMaturities(1,h+1:end).*yields(1+h:end,1)'...
                    - annualMaturities(1,1:end-h).*ExpYields_h(1:end-h,1)')/h ...
                    - repmat(yields(1,1)/12,1,numYields-h);
    xhrModel(t,:) = [nan(1,h),tmp];
    if i == 1 && modelAll.modelNum == 1
        % We use the exact expression
        xhrModel(t,:) = getExpectedExcessReturns_h1(modelAll,theta2,x,z);
    end
end

%% regressions of realized xhr on expected xhr
beta         = zeros(2,numYields);
betaRecExp   = zeros(4,numYields);
R2v          = zeros(1,numYields);
R2vRecExp    = zeros(1,numYields);
for i=h+1:numYields
    if any(isnan(xhrData(1:T-h,i)))
        startObs = find(isnan(xhrData(1:T-h,i)));
        startObs = startObs(end)+1;
    else
        startObs = 1;
    end
    % Unconditional
    Y = xhrData(startObs:T-h,i);
    X = [ones(T-h-startObs+1,1) xhrModel(startObs:T-h,i)];
    [beta(:,i),~,R2v(i),~]=robustOLS(Y,X,0,1); %h=0 lags, 1 is Newey-West se
    % Conditional on NBER recessions
    X_RecExp = [setup.NBERrec(startObs:T-h,1) setup.NBERrec(startObs:T-h,1).*xhrModel(startObs:T-h,i)...
             (1-setup.NBERrec(startObs:T-h,1)) (1-setup.NBERrec(startObs:T-h,1)).*xhrModel(startObs:T-h,i)];
    [betaRecExp(:,i),~,R2vRecExp(i),~]=robustOLS(Y,X_RecExp,0,1); %h=0 lags, 1 is Newey-West se
end


%% Output
out.xhrData      = xhrData;
out.xhrModel     = xhrModel;
out.maturities   = modelAll.matSelect;
out.beta         = beta;
out.R2v          = R2v;
out.betaRecExp   = betaRecExp;
out.R2vRecExp    = R2vRecExp;
out.numSim       = numSim;